#include "evaluate.h"


using namespace std;





// mctracks[_track][_hit] is the index of hit _hit in monte carlo track _track
//recotracks[_track][_hit] is the index of hit _hit in reconstructed track _track
// an mctrack is considered to be reconstructed if there is a recotrack with the majority of its hits coming from that mctrack, and at least required_hits hits from that track
// reco_used[_track] gets set to true if reconstructed track _track is matched to an mc track
// contributions is a map with key being number of hits from the greatest contributor, and value being the number of times such a contributor was found
void calculateEfficiency(vector<vector<unsigned int> >& mctracks, vector<vector<unsigned int> >& recotracks, unsigned int required_hits, unsigned int& nreco, vector<bool>& mc_is_reconstructed, vector<bool>& reco_used, vector<float>& mc_mom, vector<float>& reco_mom, TH1D& mom_res, map<unsigned int, unsigned int>& contributions)
{
  map<unsigned int, unsigned int> mc_hit_track;//hit index, mctrack # (# is just position in mctracks vector)
  for(unsigned int m=0;m<mctracks.size();m++)
  {
    for(unsigned int h=0;h<mctracks[m].size();h++)
    {
      mc_hit_track[mctracks[m][h]] = m;
    }
  }
  
  mc_is_reconstructed.clear();
  mc_is_reconstructed.assign(mctracks.size(), false);
  
  
  for(unsigned int r=0;r<recotracks.size();r++)
  {
    map<unsigned int, unsigned int>::iterator it;
    vector<unsigned int> mcids;//for each hit, the monte carlo track which contribute to this reco track
    for(unsigned int h=0;h<recotracks[r].size();h++)
    {
      it = mc_hit_track.find(recotracks[r][h]);
      if(it != mc_hit_track.end())
      {
        mcids.push_back(it->second);
      }
    }
    if(mcids.size() < 2){continue;}
    
    map<unsigned int, unsigned int> mclist;//mctrack #, number of hits from that track in this reco track
    map<unsigned int, unsigned int>::iterator it2;
    for(unsigned int m=0;m<mcids.size();m++)
    {
      it2 = mclist.find(mcids[m]);
      if(it2 != mclist.end())
      {
        mclist[mcids[m]]++;
      }
      else
      {
        mclist[mcids[m]] = 1;
      }
    }
    
    //find the greatest contributing mc track to the current reco track
    it2 = mclist.begin();
    unsigned int greatest = it2->first;
    unsigned int nhits = it2->second;
    for(it2 = mclist.begin();it2!=mclist.end();it2++)
    {
      if(it2->second > nhits)
      {
        greatest = it2->first;
        nhits = it2->second;
      }
    }
    if(nhits >= required_hits)
    {
      mc_is_reconstructed[greatest] = true;
      reco_used[r] = true;
      
      mom_res.Fill( (reco_mom[r] - mc_mom[greatest])/(mc_mom[greatest]) );
      
      map<unsigned int, unsigned int>::iterator mit = contributions.find(nhits);
      if(mit == contributions.end()){contributions[nhits] = 0;}
      else{contributions[nhits] += 1;}
    }
  }
  
  
  nreco = 0;
  for(unsigned int m=0;m<mc_is_reconstructed.size();m++)
  {
    if(mc_is_reconstructed[m] == true)
    {
      nreco++;
    }
  }
}


